function [top1Ts, top3Ts, objTr, objTs,time, iter] = LEML( X_tr, Y_tr, X_ts, Y_ts,k, maxIter )
%LEML Summary of this function goes here
%   Detailed explanation goes here
% obj: min||Y - XZ||^2 + 1/2(||W||^2+||H||^2)
%  Y: N*L, X: N*D , Z: D*L.  Z = WH'

%%cordinate accent
%%% estimated rank
lambda = 0;
[N,L]  =size(Y_tr);
D = size(X_tr,2);
W = randn(D,k);
w = W(:);
H = randn(L,k);
% maxIter  =20;
maxiter  = 500;
top1Ts = zeros(maxIter,1);
top3Ts = zeros(maxIter,1);
objTr = zeros(maxIter,1);
objTs = zeros(maxIter,1);
time = zeros(maxIter,1);


 
for iter = 1: maxIter
   tstart = tic;
    %%update H---closed form solution 1/2||Y - XWH'||^2+1/2 lambda*||H||^2
        A = X_tr*W;
        H = pinv(A'*A + lambda*eye(k))*A'*Y_tr;
        H = H';
    %%update w
%     W = zeros(D,k);
%     w = W(:);
    g = computegrad(X_tr,Y_tr,H,W,w,lambda);
    r=-g;d=r;
    rsold=r'*r;
    
    for i=1:maxiter
         Hd = computeHd(X_tr,H,lambda,d, D, k);%%%%%%%%%%%%%%%%%%HV
         alpha=rsold/(d'*Hd);
         w=w+alpha*d;
         r=r-alpha*Hd;
         rsnew=r'*r;
         rmse = sqrt(rsnew/size(r,1));
         if rmse<1e-4
              break;
         end
         d=r+rsnew/rsold*d;
         rsold=rsnew; 
         W = reshape(w,[D,k]);
    end
    display(i);
    temp(iter) = toc(tstart);
    W_est = W*H';
   %%%Evaluation on the training data%%%%%%%%%%%
    %%TOP 1 RATE AND TOP3 RATE
    Y_estTr = X_tr*W_est;
    objTr(iter) = norm(Y_tr - Y_estTr, 'fro')^2;
    % [top1Tr, top3Tr] = topRate(Y_estTr, Y_tr);
    % % aucTr = AUC(Y_estTr, Y_tr);
    if iter>1
        if abs(objTr(iter)-objTr(iter-1))<1e-3
            break;
        end
    end
        
    
    %%%Evaluation on the test data%%%%%%%%%%%
    %%TOP 1 RATE AND TOP3 RATE
    Y_estTs = X_ts*W_est;
    [top1Ts(iter), top3Ts(iter)] = topRate(Y_estTs, Y_ts);
    objTs(iter) = sum((Y_ts(:)-Y_estTs(:)).^2);
    if iter>3
        if objTs(iter)>objTs(iter-1) && objTs(iter-1)>objTs(iter-2) && objTs(iter-2)>objTs(iter-3)
            break;
        end
    end
    
    %  aucTs = AUC(Y_estTs, Y_ts);
end
   
for i = 1:iter
    time(i) = sum(temp(1:i));
end



end

%%%compute gradient 
function [g] = computegrad(X,Y,H,W,w,lambda)
A = X*W;
B = Y*H;
M = H'*H;
temp = X'*(A*M - B);
g = temp(:)+lambda*w;
end

%%%compute Hessian matrix
function [ Hs ] = computeHd(X,H,lambda,s, D, k)
%COMPUTEHD Summary of this function goes here
%   Detailed explanation goes here
S = reshape(s,[D,k]);
A = X*S;
M = H'*H;
temp = X'*(A*M);
Hs = temp(:)+lambda*s;

end



 

